Beads community data

Author

Masatoshi Katabuchi

Published

November 8, 2024

trait <- read_csv("data/dummy_trait.csv")
Rows: 77 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): sp
dbl (8): lma, ll, a_mass, r_mass, n_mass, p_mass, wd, sm

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
tree <- read.tree("data/dummy_tree.newick")

1 Dummy data

1.1 Habitat

x <- seq(1, 6)
y <- seq(1, 10)
z <- matrix(numeric(length(x) * length(y)), ncol = length(y))
z[1,8:9] <- z[2,8:9] <- z[3,8:9] <- 1
z[4,2:4] <- z[5,2:4] <- z[6,2:4] <- 1

z2 <- matrix(numeric(60), ncol = length(y))
z2[,1:5] <- 1

env <- tibble(env1 = as.numeric(z),
              env2 = as.numeric(z2)) %>%
  mutate(row = rep(1:6, 10)) %>%
  mutate(col = rep(1:10, each = 6)) %>%
  mutate(site = str_c("site", 1:60)) %>%
  mutate(topo = ifelse(env1 == 1, "Ridge", "valley")) %>%
  mutate(type = ifelse(env2 == 1, "Rainforest", "Rubber"))

ggplot(env, aes(y = row - 0.5, x = col - 0.5, fill = topo)) +
    geom_raster(alpha = 0.8)  +
    geom_vline(xintercept = 0:10) +
    geom_hline(yintercept = 0:6)

ggplot(env, aes(y = row - 0.5, x = col - 0.5, fill = type)) +
    geom_raster(alpha = 0.8)  +
    geom_vline(xintercept = 0:10) +
    geom_hline(yintercept = 0:6)

n_sp <- nrow(trait)
#beta_LMA <- rnorm(n_sp, log(LMA), 0.2)
beta_LMA <- 1.5
beta_WD <- 3

ts <- trait[,-1] %>% log %>% scale %>% as_tibble

1.2 Continuous environmental variables

  • soil N < 0.04%
  • soil P < 0.02%
  • moisture: moisture index (unitless)
    • large value indicates high moisture
  • spatial autocorrelation
  • but not important for species distribution

11 * 7

generate_env <- function(seed = 123, model = c("ridge_negative", "ridge_positive", "random")) {
  set.seed(seed)
  x <- seq(1, 13)
  y <- seq(1, 9)
  z <- matrix(rnorm(length(x) * length(y)), ncol = length(y))
  if (model == "ridge_positive") {
    z[3:6, 6:9] <- 2
    z[9:11, 1:4] <- 2
  } else if (model == "ridge_negative") {
    z[3:6, 6:9] <- -5
    z[9:11, 1:4] <- -5
  }

  for (i in 2:(length(x)-1)) {
    for (j in 2:(length(y)-1)) {
      z[i, j] <-
        rnorm(1,
              mean(
                c(z[i - 1, j - 1],
                z[i - 1, j],
                z[i - 1, j + 1],
                z[i, j - 1],
                z[i, j + 1],
                z[i + 1, j - 1],
                z[i + 1, j],
                z[i + 1, j + 1])),
              0)
    }
  }
  z <- z[,-9]
  z <- z[,-1]
  z <- z[-13,]
  z <- z[-1,]

  x <- seq(1, 11)
  y <- seq(1, 7)
  x2 <- rep(x, length(y))
  y2 <- rep(y, each = length(x))
  # as.numeric(z)
  # as.numeric(t(z))

  dat <- tibble(x = x2 |> as.numeric(),
                y = y2,
                z = as.numeric(z) )
    # mutate(z = ifelse(x > 5, z + 1, z))

  dat2 <- dat %>%
    mutate(x = x - 1) %>%
    mutate(y = y - 1)
  dat2
}

moist <- generate_env(model = "ridge_negative")
nitro <- generate_env(model = "random")
phos <- generate_env(model = "random")

nitro <- nitro |>
  mutate(soil_n = exp(z)) |>
  mutate(soil_n = soil_n / max(soil_n) * 0.4)

phos <- phos |>
  mutate(soil_p = exp(z)) |>
  mutate(soil_p = soil_p / max(soil_p) * 0.2)

moist <- moist |>
  mutate(soil_moist = z)

dat2 <- nitro |>
  dplyr::select(-z) |>
  mutate(soil_p = phos$soil_p) |>
  mutate(soil_moist = moist$soil_moist)

ggplot(dat2, aes(x = x, y = y, fill = soil_n |> log())) +
  geom_raster()

ggplot(dat2, aes(x = x, y = y, fill = soil_p |> log())) +
  geom_raster()

ggplot(dat2, aes(x = x, y = y, fill = soil_moist)) +
  geom_raster()

tmp_fun <- function(data, env) {
  tmp <- NULL
  for (i in 0:9) {
    for (j in 0:5) {
       s1 <- data |>
          filter(x == i & y == j) |>
          pull({{env}})
       s2 <- data |>
          filter(x == (i + 1) & y == j) |>
          pull({{env}})
       s3 <- data |>
          filter(x == i  & y == (j + 1)) |>
          pull({{env}})
       s4 <- data |>
          filter(x == i  & y == (j + 1)) |>
          pull({{env}})
      tmp <- c(tmp, mean(c(s1, s2, s3, s4)))
    }
  }
  tmp
}

env3 <- env |>
  mutate(soil_moist = tmp_fun(dat2, soil_moist))

ggplot(env3, aes(x = col, y = row, fill = soil_moist)) +
  geom_raster()

dat2 |>
   dplyr::select(x, y, soil_n, soil_p, soil_moist) |>
   round(3) |>
   write_csv("data/soil.csv")

1.3 Abundance

Abundance of ith species at site j:

\[ Abundance_{ij} \sim Pois(\lambda_{ij}) \] \[ \log\lambda_{ij} = \beta_{0j} + \beta_{1j} \times Ridge_i + \beta_{2j} \times Forest_i + \beta_{3j} \times Ridge_i \times Forest_i \] \[ \beta_{mj} = \gamma_{m0} + \gamma_{m1} \times logLMA_j + \gamma_{m2} \times logWD_j + \epsilon_{mj} \] \[ \epsilon_{mj} \sim N(0, \sigma_m) \]

beta0 <- rnorm(n_sp, 5, 0.2)

beta1 <- rnorm(n_sp,
               mean = -5 + 1.5 * ts$lma + 0.6 * ts$wd,
               sd = 0.2)

beta2 <- rnorm(n_sp, -1 + 0.8 * ts$sm, 0.2)

beta3 <- rnorm(n_sp,
               mean = 1 * ts$lma,
               sd = 0.2)

beta4 <- rnorm(n_sp,
               mean = -0.5 * ts$lma,
               sd = 0.2)

beta <- cbind(beta0, beta1, beta2, beta3, beta4)

env2 <- tibble(intercept  = 1,
               env1 = env$env1,
               env2 = env$env2,
               interaction = env$env2 * env$env1,
               moist = env3$soil_moist)

env2
# A tibble: 60 × 5
   intercept  env1  env2 interaction   moist
       <dbl> <dbl> <dbl>       <dbl>   <dbl>
 1         1     0     1           0  0.193 
 2         1     0     1           0  0.0693
 3         1     0     1           0 -0.280 
 4         1     0     1           0 -1.18  
 5         1     0     1           0 -1.86  
 6         1     0     1           0 -2.41  
 7         1     0     1           0  0.198 
 8         1     0     1           0  0.0913
 9         1     0     1           0 -0.652 
10         1     1     1           1 -1.94  
# ℹ 50 more rows
summary(ts)
      lma                 ll              a_mass             r_mass        
 Min.   :-2.35492   Min.   :-1.9081   Min.   :-2.63320   Min.   :-2.86039  
 1st Qu.:-0.77439   1st Qu.:-0.6445   1st Qu.:-0.51068   1st Qu.:-0.65759  
 Median :-0.07205   Median :-0.1391   Median :-0.01897   Median : 0.03553  
 Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.00000  
 3rd Qu.: 0.73571   3rd Qu.: 0.4938   3rd Qu.: 0.64283   3rd Qu.: 0.64807  
 Max.   : 2.91324   Max.   : 2.8215   Max.   : 2.39729   Max.   : 2.05776  
     n_mass            p_mass               wd                 sm          
 Min.   :-2.7071   Min.   :-2.65152   Min.   :-3.20202   Min.   :-2.14566  
 1st Qu.:-0.8276   1st Qu.:-0.69644   1st Qu.:-0.73935   1st Qu.:-0.65216  
 Median : 0.1815   Median : 0.05989   Median :-0.01349   Median :-0.02183  
 Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.00000  
 3rd Qu.: 0.7325   3rd Qu.: 0.67051   3rd Qu.: 0.62138   3rd Qu.: 0.58372  
 Max.   : 2.1057   Max.   : 1.92002   Max.   : 2.15467   Max.   : 2.51892  
mu <- beta %*% t(env2)
ab <- matrix(numeric(60*77), ncol = ncol(mu))
for (i in 1:nrow(mu)) {
  for (j in 1:ncol(mu)) {
    ab[i, j] <- rpois(1, exp(mu[i, j]))
    #ab[i, j] <- rpoilog(1, exp(mu[i, j]), 1)
  }
}

2 Species abundance

# Function to plot abundance for top/bottom species based on a specific trait
plot_abundance <- function(env, ab, species_list, title_list) {
  for (i in seq_along(species_list)) {
    p <- env %>%
      mutate(abund = ab[species_list[i], ]) %>%
      ggplot(aes(x = col - 0.5, y = row - 0.5, fill = abund)) +
      ggtitle(title_list[i]) +
      geom_raster() +
      geom_vline(xintercept = 0:10) +
      geom_hline(yintercept = 0:6)
    print(p)
  }
}

# Set rownames for abundance matrix `ab`
rownames(ab) <- trait$sp

2.1 Plot Top 5 LMA

top_sp_lma <- trait %>%
  arrange(-lma) %>%
  slice(1:5) %>%
  pull(sp)
plot_abundance(env, ab, top_sp_lma, top_sp_lma)

2.2 Plot Lowest 5 LMA

low_sp_lma <- trait %>%
  arrange(lma) %>%
  slice(1:5) %>%
  pull(sp)
plot_abundance(env, ab, low_sp_lma, low_sp_lma)

2.3 Plot Top 5 Seed Mass

top_sp_sm <- trait %>%
  arrange(-sm) %>%
  slice(1:5) %>%
  pull(sp)
plot_abundance(env, ab, top_sp_sm, top_sp_sm)

2.4 Plot Lowest 5 Seed Mass

low_sp_sm <- trait %>%
  arrange(sm) %>%
  slice(1:5) %>%
  pull(sp)

plot_abundance(env, ab, low_sp_sm, low_sp_sm)

3 Species list

3.1 Abundant species or species with clear patterns

abundant_sp <- c(top_sp_lma, low_sp_lma, top_sp_sm, low_sp_sm) |>
  unique() |>
  sort()
abundant_sp
 [1] "Beilschmiedia_robusta"    "Camellia_forrestii"      
 [3] "Camellia_Pitardii"        "Castanopsis_wattii"      
 [5] "Claoxylon_khasianum"      "Eurya_trichocarpa"       
 [7] "Ficus_henryi"             "Ficus_heteromorpha"      
 [9] "Itoa_orientalis"          "Lindera_caudata"         
[11] "Lithocarpus_crassifolius" "Lithocarpus_henryi"      
[13] "Lithocarpus_xylocarpus"   "Machilus_salicina"       
[15] "Neolitsea_polycarpa"      "Rhododendron_decorum"    
[17] "Rhus_sylvestris"          "Schima_noronhae"         
[19] "Symplocos_anomala"       
plot_abundance(env, ab, abundant_sp, abundant_sp)

3.2 Other species

rare_sp <- trait %>%
  filter(!(sp %in% abundant_sp)) %>%
  pull(sp) |>
  # setdiff(trait$sp) |>
  sort()

rare_sp
 [1] "Acer_campbellii"             "Actinodaphne_forrestii"     
 [3] "Alnus_nepalensis"            "Anneslea_fragrans"          
 [5] "Betula_alnoides"             "Camellia_assamica"          
 [7] "Camellia_yunnanensis"        "Castanopsis_echidnocarpa"   
 [9] "Craibiodendron_yunnanense"   "Cyclobalanopsis_stewardiana"
[11] "Eriobotrya_bengalensis"      "Eriobotrya_tengyuehensis"   
[13] "Eurya_paratetragonoclada"    "Eurya_quinquelocularis"     
[15] "Exbucklandia_populnea"       "Ficus_racemosa"             
[17] "Glochidion_a"                "Ilex_corallina"             
[19] "Ilex_manneiensis"            "Illicium_macranthum"        
[21] "Laurocerasus_undulata"       "Leucosceptrum_canum"        
[23] "Lindera_thomsonii"           "Lithocarpus_confinis"       
[25] "Lithocarpus_dealbatus"       "Lithocarpus_hancei"         
[27] "Lithocarpus_truncatus"       "Litsea_a"                   
[29] "Litsea_chunii"               "Litsea_elongata"            
[31] "Machilus_gamblei"            "Machilus_yunnanensis"       
[33] "Mahonia_duclouxiana"         "Manglietia_insignis"        
[35] "Melia_toosendan"             "Meliosma_arnottiana"        
[37] "Meliosma_kirkii"             "Michelia_floribunda"        
[39] "Myrsine_semiserrata"         "Neolitsea_chuii"            
[41] "Olea_rosea"                  "Osbeckia_mepalensis"        
[43] "Padus_napaulensis"           "Rhododendron_irroratum"     
[45] "Rhododendron_leptothrium"    "Saurauia_napaulensis"       
[47] "Schefflera_fengii"           "Schefflera_shweliensis"     
[49] "Schima_argentea"             "Skimmia_arborescens"        
[51] "Sterculia_nobilis"           "Stewartia_pteropetiolata"   
[53] "Stewartia_yunnanensis"       "Symplocos_grandis"          
[55] "Symplocos_ramosissima"       "Symplocos_sumuntia"         
[57] "Ternstroemia_gymnanthera"    "Vaccinium_duclouxii"        
tibble(sp = abundant_sp, abund = "abundant") |>
  bind_rows(tibble(sp = rare_sp, abund = "rare")) |>
  mutate(sp2 = str_replace_all(sp, "_", " ")) |>
  write_csv("data/note.csv")
plot_abundance(env, ab, rare_sp, rare_sp)